Unsupervised Lab

Data Description

gdat is Gene Expression Data, n = 445 patients x p = 353 genes

• Only 353 genes with somatic mutations from COSMIC are retained ## Data is Level III TCGA BRCA RNA-Sequencing gene expression data that have already been pre-processed according to the following steps:
• Reads normalized by RPKM
• Corrected for overdispersion by a log-transformation (1 + data)
• Short gene name labels are given as the column names

cdat is Clinical Data, n = 445 patients x q = 6 clinical features

• Subtype - denotes 5 PAM50 subtypes including Basal-like, Luminal A, Luminal B, HER2-enriched, and Normal-like
• ER-Status - estrogen-receptor status
• PR-Status - progesterone-receptor status
• HER2-Status - human epidermal growth factor receptor 2 status
• Node - number of lymph nodes involved
• Metastasis - indicator for whether the cancer has metastasized

Problems

Problem 1 - Dimension reduction

1a - Apply PCA, NMF, ICA and MDS, UMAP, and tSNE to this dataset. Compare and contrast the results using these methods.
1b - Relate the dimension reduction results with the clinical data. Is any clinical information reflected in the lower dimensional spaces?

1c - Overall, which dimension reduction method do you recommend for this data set and why?

Problem 2 - Clustering

2a - Apply various clustering algorithms such as K-means (explore different K), hierarchical clustering (explore different linkages), NMF, and biclustering. Compare the clustering results using these methods.

2b - Relate the clustering results with the clinical data. Can the clustering algorithm recover some of the clinical information such as cancer subtypes?

2c - Use Consensus Clustering to help Validate Clustering Results

2c - Overall, which clustering method(s) do you recommend for this data set and why?

Problem 3 - Multiple comparisons

3a - Identify important genes to differetiate ER postive and negative, PR postive and negative, HER2 postive and negative, metastasis status.

3b - Try different procedures to adjust for multiple comparisons.

3c - Examine the lists of genes identified using different methods for each clinical response. Which method is best? Why?

Problem 5 - Graphical models

5a - Use graphical models to explore interactions among genes. Are any of the well-connected genes related to patterns previously identified?

Problem 6 - Visulaization

6a - Visualize this data using multiple approaches.

6b - Prepare the “best” visual summary of this data.

Problem 7 - Exploratory Data Analysis Summary

7a - What is the most interesting finding?

7b - Is this finding consistent and stable?

7c - Prepare a visual summary that best illustrates this interesting finding.

R scripts to help out with the BRCA case study Lab Don’t peek at this if you want to practice coding on your own!!

Load Data

load("UnsupL_SISBID_2020.Rdata")
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.6.2

library(ConsensusClusterPlus)
library(kknn)
library(GGally)

## Warning: package 'GGally' was built under R version 3.6.2

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library(umap)

## Warning: package 'umap' was built under R version 3.6.2

library(Rtsne)
library(igraph)

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

library(XMRF)
library(spacejam)

## Loading required package: splines

## Loading required package: Matrix

Explore Data

dim(gdat)

## [1] 445 353

dim(cdat)

## [1] 445   6

table(cdat$Subtype)

## 
##    Basal-like HER2-enriched     Luminal A     Luminal B   Normal-like 
##            79            53           200           106             7

table(cdat$ER)

## 
##               Indeterminate                    Negative 
##                           2                         100 
##               Not Performed Performed but Not Available 
##                           2                           2 
##                    Positive 
##                         339

table(cdat$PR)

## 
##               Indeterminate                    Negative 
##                           3                         147 
##               Not Performed Performed but Not Available 
##                           2                           2 
##                    Positive 
##                         291

table(cdat$HER2)

## 
##     Equivocal      Negative Not Available      Positive 
##             5           370             4            65

table(cdat$Node)

## 
##   0   1   2   3 
## 221 146  54  23

table(cdat$Metastasis)

## 
##   0   1 
## 427  11

table(cdat$ER,cdat$PR)

##                              
##                               Indeterminate Negative Not Performed
##   Indeterminate                           0        1             0
##   Negative                                1       93             0
##   Not Performed                           0        0             2
##   Performed but Not Available             0        0             0
##   Positive                                2       53             0
##                              
##                               Performed but Not Available Positive
##   Indeterminate                                         0        1
##   Negative                                              0        6
##   Not Performed                                         0        0
##   Performed but Not Available                           2        0
##   Positive                                              0      284

#cluster heatmap - biclustering

#cluster heatmap - biclustering
aa = grep("grey",colors())
bb = grep("green",colors())
cc = grep("red",colors())
gcol2 = colors()[c(aa[1:2],bb[1:25],cc[1:50])]

Without scaling

heatmap(gdat,col=gcol2,hclustfun=function(x)hclust(x,method="ward.D"))

With scaling

heatmap(scale(gdat),col=gcol2,hclustfun=function(x)hclust(x,method="ward.D"))

Cols=function(vec){cols=rainbow(length(unique(vec)))
                   return(cols[as.numeric(as.factor(vec))])}

heatmap(scale(gdat),col=gcol2,hclustfun=function(x)hclust(x,method="ward.D"),labRow=cdat$Subtype,RowSideColors=Cols(cdat$Subtype))

#Dimension Reduction

PCA

Cols=function(vec){cols=rainbow(length(unique(vec)))
                   return(cols[as.numeric(as.factor(vec))])}

sv = svd(scale(gdat,center=TRUE,scale=FALSE))
V = sv$v
Z = gdat%*%V

K = 3
pclabs = c("PC1","PC2","PC3","PC4")
par(mfrow=c(1,K))
for(i in 1:K){
  j = i+1
  plot(Z[,i],Z[,j],pch=16,xlab=pclabs[i],ylab=pclabs[j],col=Cols(cdat$Subtype))
}
legend(-45,0,pch=16,col=rainbow(5),levels(cdat$Subtype))

Pairs Plot

PC1<-as.matrix(Z[,1])
PC2<-as.matrix(Z[,2])
PC3<-as.matrix(Z[,3])
PC4<-as.matrix(Z[,4])
PC5<-as.matrix(Z[,5])

pc.df.cdat<-data.frame(Subtype = cdat$Subtype, PC1, PC2, PC3, PC4, PC5)
ggpairs(pc.df.cdat, mapping = aes(color = Subtype))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

MDS

Dmat = dist(gdat,method="maximum")
mdsres = cmdscale(Dmat,k=2)
plot(mdsres[,1],mdsres[,2],pch=16,col=Cols(cdat$Subtype), main = "Dimension Reduction MDS")
legend(-30,20,pch=16,col=rainbow(5),levels(cdat$Subtype))

ICA

require("fastICA")

## Loading required package: fastICA

K = 4
icafit = fastICA(gdat,n.comp=K)

kk = 3
pclabs = c("ICA1","ICA2","ICA3","ICA4")
par(mfrow=c(1,kk))
for(i in 1:kk){
  j = i+1
  plot(icafit$A[i,],icafit$A[j,],pch=16,xlab=pclabs[i],ylab=pclabs[j],col=Cols(cdat$Subtype))
}
legend(-1,2.8,pch=16,col=rainbow(5),levels(cdat$Subtype))

UMAP

gdat.umap = umap(gdat)
plot(gdat.umap$layout[,1], y =gdat.umap$layout[,2], type = "n", main = "UMAP", xlab = "UMAP1", ylab = "UMAP2")
text(gdat.umap$layout[,1], y =gdat.umap$layout[,2], type = "n", cdat$Subtype, col=Cols(cdat$Subtype), cex = .7 )

## Warning in text.default(gdat.umap$layout[, 1], y = gdat.umap$layout[, 2], :
## graphical parameter "type" is obsolete

#Clustering

K-means

K = 5
km = kmeans(gdat,centers=K,nstart=25)
table(km$cluster,cdat$Subtype)

##    
##     Basal-like HER2-enriched Luminal A Luminal B Normal-like
##   1          5             8        45        12           4
##   2          0             2        42        60           1
##   3         74             5         1         1           1
##   4          0             7       109        26           1
##   5          0            31         3         7           0

Plot Kmeans with labels

plot(Z[,1],Z[,2],col=km$cluster, main = "Plot Kmeans Clusters ", xlab = "PC1", ylab = "PC2")
text(Z[,1],Z[,2],cdat$Subtype,cex=.75,col=km$cluster)
cens = km$centers
points(cens%*%V[,1],cens%*%V[,2],col=1:K,pch=16,cex=3)

Hierarchical

#which linakge is the best?
#which distance metric is the best?

Dmat = dist(gdat)
com.hc = hclust(Dmat,method="ward.D")


plot(com.hc,labels=cdat$Subtype,cex=.5)

res.com = cutree(com.hc,5)
table(res.com,cdat$Subtype)

##        
## res.com Basal-like HER2-enriched Luminal A Luminal B Normal-like
##       1          1             3        95        11           3
##       2          0             4        73        65           1
##       3         75             4         5         4           1
##       4          0            27         3         7           0
##       5          3            15        24        19           2

Consensus Clustering with Hierarchical

results = ConsensusClusterPlus(gdat,maxK=6,reps=500,pItem=0.8,pFeature=1,
clusterAlg="hc",distance="pearson",plot="png")

## end fraction

## clustered
## clustered
## clustered
## clustered
## clustered

Look at results for first 5 clusters

 heatmap(results[[2]][["consensusMatrix"]][1:5,1:5])

Spectral Clustering

K = 5
s_gdat = specClust(gdat, centers=K, nn = 7, method = "symmetric", gmax=NULL)

Visualize

plot(Z[,1],Z[,2],col=s_gdat$cluster, main = "Visualize Spectral Clusters", xlab = "PC1", ylab = "PC2")
text(Z[,1],Z[,2],cdat$Subtype,cex=.75,col=s_gdat$cluster)

Genes significantly associated with ER or PR Status, etc

x = gdat[cdat$ER=="Positive" | cdat$ER=="Negative",]
y.er = cdat$ER[cdat$ER=="Positive" | cdat$ER=="Negative"]
y.label = rep(1, length(y.er))
y.label[y.er == "Positive"]=2

ps = NULL
for(i in 1:ncol(gdat)) ps = c(ps,
 t.test(x[y.label==1,i],x[y.label==2,i])$p.value)
fdrs.bh = p.adjust(ps, method="BH")

cat("Number of Tests significant with alpha=0.1 using Bonferroni correction:",
sum(ps<0.1/length(y.label)), fill=TRUE)

## Number of Tests significant with alpha=0.1 using Bonferroni correction: 165

cat("Number of Tests with FDR below 0.1:",
sum(fdrs.bh<0.1), fill=TRUE)

## Number of Tests with FDR below 0.1: 259

plot(sort(ps,decreasing=FALSE),ylab="P-Values")
#BH procedure
abline(a=0, b=0.1/length(y.label),col="red")
#Bonferroni
abline(a=0.1/length(y.label), b=0,col="blue")

Graphical models - how are genes related?

#this takes a bit of time

lam = lambdaMax(gdat)*sqrt(log(ncol(gdat))/nrow(gdat))*0.01
net = XMRF(t(gdat),method="LPGM",lambda.path=lam,N=1,th=.0025)

## LPGM ::: Conducting sampling ... in progress:  100 % 
## LPGM Completed.

plot(net)

## Plot optimal network

##                [,1]     [,2]
##   [1,] 19.838264657 30.86689
##   [2,]  8.295297756 30.98339
##   [3,] -4.945591598 48.50826
##   [4,] -0.072424102 38.72952
##   [5,]  1.446459793 35.98709
##   [6,]  3.194760332 39.58998
##   [7,] 21.312783923 30.77726
##   [8,]  8.373994622 35.18537
##   [9,]  4.707045609 57.09021
##  [10,]  5.158871953 33.55555
##  [11,] -2.730617466 37.54865
##  [12,]  5.678206097 56.22942
##  [13,] 25.438684620 40.82255
##  [14,]  8.291445725 55.49436
##  [15,]  3.042894993 56.59438
##  [16,] -1.263971134 50.76479
##  [17,] 26.142476160 45.51535
##  [18,] 12.425998424 31.38349
##  [19,]  8.502924125 41.45538
##  [20,] -0.200305414 33.08569
##  [21,] 16.081662098 51.65025
##  [22,] 21.917792208 51.52781
##  [23,] 21.832471278 47.67047
##  [24,] 26.599573683 41.47744
##  [25,]  3.565047891 37.93286
##  [26,]  4.339860804 31.75050
##  [27,]  7.324784847 52.33422
##  [28,] 20.228125422 34.05805
##  [29,]  4.852517564 31.08728
##  [30,]  3.292012143 31.96992
##  [31,] -1.607282823 40.33197
##  [32,] 11.220051034 35.32044
##  [33,] -0.480258757 53.37978
##  [34,] -1.903283851 37.86544
##  [35,] 20.999739217 43.35891
##  [36,] 23.182693792 33.01077
##  [37,]  6.863407060 57.16108
##  [38,]  9.367402156 44.12062
##  [39,]  9.377171601 49.81023
##  [40,] 18.671656093 56.29594
##  [41,]  2.724928650 34.16188
##  [42,]  7.612850869 41.95366
##  [43,] 18.230046824 46.32492
##  [44,]  0.660423176 38.95581
##  [45,] 13.898747569 35.62331
##  [46,] 13.841358863 56.75667
##  [47,] 19.540657234 47.48261
##  [48,]  1.343899473 32.35624
##  [49,] 10.834075224 56.98855
##  [50,] 22.109606361 53.52758
##  [51,] -3.960265652 47.45734
##  [52,]  1.569623673 56.45270
##  [53,] 25.932057403 37.54140
##  [54,]  1.544511485 37.29487
##  [55,] 10.765269363 34.83927
##  [56,] 20.225106604 48.75143
##  [57,] 16.589417314 39.41186
##  [58,] 11.102320523 36.93196
##  [59,] 14.057238179 36.67035
##  [60,] -3.445348623 51.58113
##  [61,]  1.444924244 40.59768
##  [62,] 10.489475429 32.71512
##  [63,]  1.030163536 40.40300
##  [64,]  9.296411878 34.64020
##  [65,] 14.347754211 38.13310
##  [66,] 22.678879007 40.97109
##  [67,] 20.026789236 55.50454
##  [68,] 13.034252033 55.40708
##  [69,] 17.519550403 53.03269
##  [70,] 19.483331214 32.10677
##  [71,] 26.363384270 44.13467
##  [72,] 13.419331348 36.90529
##  [73,]  0.670777338 39.97927
##  [74,] -3.683297526 49.08652
##  [75,]  7.900145849 34.66698
##  [76,]  2.802066170 38.97391
##  [77,] -2.031177041 53.58478
##  [78,]  4.737252968 36.27136
##  [79,]  5.820065189 36.89361
##  [80,]  0.904722613 31.75210
##  [81,] -2.470066043 48.47391
##  [82,]  5.853995445 50.76320
##  [83,] 10.120748167 41.99061
##  [84,] 13.868630933 58.21228
##  [85,] 18.428870221 43.32143
##  [86,]  1.007264220 52.35334
##  [87,] 20.355687570 37.99708
##  [88,] -0.196485013 42.40738
##  [89,] 22.865013156 44.52788
##  [90,]  4.295255982 40.56307
##  [91,]  8.838605013 51.75406
##  [92,]  0.798961138 54.12138
##  [93,] 20.549367075 35.59592
##  [94,] 11.364135124 30.94853
##  [95,]  8.146185457 37.44643
##  [96,]  5.765998187 39.80018
##  [97,] 15.003234626 46.87180
##  [98,]  7.534188860 50.54428
##  [99,]  9.401946764 56.36881
## [100,] 10.515262743 35.83292
## [101,] 15.671676728 27.14843
## [102,] 11.398655079 55.64188
## [103,] 12.744248299 53.83570
## [104,]  1.792001410 38.67687
## [105,] 11.792116145 31.51406
## [106,]  7.784714063 31.11862
## [107,] 16.152423957 54.74328
## [108,] 21.971439597 32.44607
## [109,]  7.150080834 36.21841
## [110,]  2.463621634 51.29772
## [111,] 17.347555735 55.97895
## [112,] -1.971338480 36.86097
## [113,] 12.549164880 39.63456
## [114,] 11.156177403 58.63994
## [115,] 19.793871747 45.91049
## [116,] 14.290558433 37.73597
## [117,]  5.197463501 36.48672
## [118,] 24.671787425 50.02539
## [119,] -1.024601796 41.11683
## [120,]  1.546314270 39.72496
## [121,] 15.240984587 57.90895
## [122,] 17.516320397 50.84673
## [123,] 13.262042582 39.62001
## [124,]  3.558381286 36.67464
## [125,]  3.088762484 33.63772
## [126,] 13.751097631 51.41151
## [127,]  9.807747524 54.80725
## [128,]  4.658885303 44.19630
## [129,] 11.029461087 33.43757
## [130,]  4.504082787 33.56112
## [131,]  3.689740068 36.99346
## [132,]  4.016149566 42.37604
## [133,] 24.062312752 40.30369
## [134,]  7.282962721 58.66570
## [135,]  5.542570977 41.10222
## [136,]  5.164241054 31.68914
## [137,] 21.361638072 44.86243
## [138,] 10.443698628 52.73173
## [139,]  4.617545296 39.94560
## [140,]  1.805002569 41.21817
## [141,] -0.625038606 43.38235
## [142,]  5.922463717 53.17410
## [143,] -0.795226171 33.36801
## [144,]  9.505580830 36.05210
## [145,]  5.947356835 41.23146
## [146,]  6.741400740 43.88794
## [147,]  6.232160576 38.06550
## [148,] 18.174144724 29.47499
## [149,] 18.984864842 29.75727
## [150,]  8.194926084 57.56744
## [151,]  2.737755542 37.92951
## [152,]  3.566641382 40.61698
## [153,] 13.421756629 35.27450
## [154,] 12.961315673 38.91233
## [155,]  2.201081793 34.97813
## [156,] 14.300970591 34.29035
## [157,] 14.522505545 35.91764
## [158,] 14.194409995 35.07020
## [159,]  5.794803148 33.80872
## [160,] 23.048211645 39.16700
## [161,] 13.221070178 36.27409
## [162,]  2.173762497 53.82266
## [163,]  7.131777364 33.46496
## [164,] 22.602932716 31.35726
## [165,] -1.242489866 42.85625
## [166,] 24.184834115 38.28010
## [167,] 11.058047615 30.40826
## [168,] 16.914687337 28.00923
## [169,]  0.595773901 36.44548
## [170,]  8.727065121 36.19884
## [171,]  3.503497471 31.44125
## [172,] 23.008958362 52.38288
## [173,]  3.502070291 36.37646
## [174,]  1.088052448 38.82101
## [175,] 13.852374093 36.18074
## [176,] 12.271326498 35.70071
## [177,] 11.831185597 33.42746
## [178,]  6.627087544 34.57717
## [179,] -5.305666592 46.59879
## [180,] 13.528831284 47.98797
## [181,] 24.605389472 46.43136
## [182,] 17.020539887 57.39635
## [183,]  9.426264208 31.30048
## [184,]  1.355908408 35.27176
## [185,] 12.231785935 52.32198
## [186,]  9.420958674 38.06781
## [187,]  8.943182968 53.49179
## [188,] 26.394092588 39.62829
## [189,] 14.829038658 55.71897
## [190,]  6.832353517 55.45470
## [191,] -2.358970474 39.54544
## [192,] -1.666506659 34.04969
## [193,] -2.275242197 52.24315
## [194,]  0.339304584 50.77731
## [195,]  8.551931395 35.66876
## [196,] 21.490377333 39.61240
## [197,] 25.301002346 48.27629
## [198,]  4.008223791 50.72216
## [199,] 25.306556265 38.95048
## [200,]  8.860426188 30.77092
## [201,] -0.281277037 35.18854
## [202,]  0.539689078 35.22785
## [203,] 14.283730116 52.82004
## [204,] 16.446448589 46.19847
## [205,] 18.509402870 51.84716
## [206,] -1.287485890 36.21155
## [207,]  0.532260756 40.77345
## [208,]  8.780027709 31.58064
## [209,]  2.132271940 46.57556
## [210,]  2.905911066 44.07795
## [211,] -0.117219586 32.40754
## [212,]  9.109844984 37.12172
## [213,]  4.884630741 34.96444
## [214,]  3.990345452 35.60359
## [215,]  6.306044925 39.00490
## [216,] 20.121503415 51.39196
## [217,] 21.567968334 36.46362
## [218,]  9.346722849 33.57526
## [219,]  8.480829085 39.30396
## [220,]  1.805003731 37.88845
## [221,]  8.942536751 40.12915
## [222,] 12.327058669 57.07843
## [223,]  0.530713842 37.01683
## [224,] -0.202187719 37.73361
## [225,]  5.881059705 39.26675
## [226,] -0.434982039 36.91457
## [227,] 20.754545431 32.76185
## [228,] -2.549200370 50.29915
## [229,]  7.427969542 37.70748
## [230,] 13.794024250 39.06348
## [231,] 26.163350380 42.71873
## [232,] 20.882350744 52.73482
## [233,]  6.298537381 42.66674
## [234,]  3.075718331 35.74845
## [235,]  9.117956833 40.45187
## [236,] 22.289352099 42.49072
## [237,]  9.378959486 59.03277
## [238,]  0.461578119 32.35721
## [239,] 11.962596039 50.65777
## [240,]  1.629106732 36.40548
## [241,] 15.930469830 56.58868
## [242,] 17.370782612 44.81138
## [243,]  9.741241772 31.83634
## [244,] 15.422347606 32.64791
## [245,]  2.992743180 52.61110
## [246,] 13.324546472 49.61528
## [247,]  0.861906584 37.36720
## [248,]  3.746359467 35.07027
## [249,]  4.883874289 52.09552
## [250,]  3.829502191 33.80065
## [251,] 23.706780285 51.17331
## [252,] 10.578815917 34.12663
## [253,]  9.855902993 32.74570
## [254,] 10.398643425 30.19314
## [255,] 24.095915888 48.75595
## [256,] 25.840009426 46.96213
## [257,] -0.672700369 52.08543
## [258,]  5.450968511 33.94353
## [259,] 21.897124164 37.98763
## [260,]  9.961617224 35.69469
## [261,] 17.661452824 54.53311
## [262,] 14.504123072 36.45517
## [263,] 23.090094348 46.01576
## [264,]  3.196178189 42.98282
## [265,]  5.678252169 31.34060
## [266,] 23.421533470 47.53541
## [267,] -0.742544464 39.61887
## [268,] 21.102362475 50.20402
## [269,] -3.429137838 41.17751
## [270,] 13.320598350 34.63415
## [271,]  7.405313263 37.45012
## [272,]  2.705068806 39.73527
## [273,] 20.442421975 29.44126
## [274,] 15.968419945 53.14040
## [275,] 11.230467717 54.14501
## [276,]  3.942458003 55.26329
## [277,] -0.116401910 39.28005
## [278,]  2.783388445 36.34478
## [279,]  0.482561786 55.56206
## [280,] 19.563960474 42.13002
## [281,]  2.520830116 55.20942
## [282,]  2.590745759 35.24679
## [283,] -0.576073795 40.37840
## [284,]  9.797458107 40.55278
## [285,]  7.948746557 41.28381
## [286,]  7.765386982 40.88162
## [287,]  5.397004706 40.58520
## [288,]  4.986185106 33.19470
## [289,]  0.760662853 38.42472
## [290,] 24.655935356 36.79929
## [291,] 24.571517037 44.93621
## [292,]  3.076881019 41.44441
## [293,] 24.923961817 43.50148
## [294,]  2.505749609 36.85567
## [295,]  4.674507267 46.79393
## [296,] 21.028445051 41.44502
## [297,] 24.247606757 33.72549
## [298,]  1.542356920 50.18198
## [299,]  7.081825746 35.44002
## [300,] 17.159343282 47.79718
## [301,] 21.936724774 34.23455
## [302,]  9.661174753 57.83958
## [303,] 19.103697846 54.49899
## [304,] 14.416119074 54.26214
## [305,] -0.909443890 49.33533
## [306,] 18.400181200 28.08469
## [307,] 15.259638461 48.61540
## [308,]  5.279185558 35.66529
## [309,]  7.389919538 32.94221
## [310,] 19.362738935 53.09857
## [311,] 23.555292973 43.15174
## [312,] 23.949805001 35.01009
## [313,] 15.980224429 33.49990
## [314,] 19.268796328 50.07388
## [315,]  2.039220549 41.85566
## [316,] -0.786454720 54.79682
## [317,] 19.534304324 44.40725
## [318,] 16.622968016 49.59964
## [319,] 22.936077684 50.00550
## [320,] 11.523671387 49.00336
## [321,] -4.301039961 50.28031
## [322,] 21.355132041 46.40529
## [323,]  4.092743574 53.56517
## [324,]  3.354968173 38.87101
## [325,]  3.629821483 57.72785
## [326,]  8.759446640 34.10995
## [327,] 12.830559872 34.96710
## [328,] 13.550685049 38.41369
## [329,] -0.005711561 43.81178
## [330,] 10.431439722 51.06879
## [331,]  5.470333110 38.35883
## [332,]  7.756134211 38.58301
## [333,]  2.228544534 35.61018
## [334,] 25.191868878 35.58392
## [335,]  7.377755388 54.05618
## [336,] 20.873025165 54.34941
## [337,]  4.107299070 31.11167
## [338,] 24.403930085 41.89686
## [339,] 14.995256574 50.43140
## [340,] 12.514656907 58.61875
## [341,]  5.681890937 58.22659
## [342,] 22.714005711 35.44925
## [343,] 23.138558363 37.01331
## [344,] 20.030065198 40.11837
## [345,]  1.231276141 43.27957
## [346,]  4.466884527 38.17857
## [347,]  7.450070180 32.09569
## [348,]  1.960576018 39.02725
## [349,]  7.389100441 43.11813
## [350,] 18.299190485 48.85301
## [351,] 22.034390427 49.00946
## [352,]  9.916337207 38.30725
## [353,]  5.277049733 54.67881

## estimation using spacejam
net2 = SJ(gdat)
net2 = net2$G
sjnet <- 1*(net2)[,,1]
sjnet = graph.adjacency(sjnet, mode="undirected")
plot(sjnet, vertex.size=0.1, vertex.label=NA)